function [dPsieffBydLamba] = PsieffDerivativeNew(lambda_i, fVec, Params)
% 
% PSIEFFDERIVATE - compute the derivative of the effective phase 
%     Psi_eff = Psi(f; M, eta, tc) 
% with respect to the re-parametrized set of parameters {M, eta, A, B, a1, a2, a3, a4, tc}
% usage: [dPsieffBydLamba] = PsieffDerivativeNew(lambda_i, fVec, Params)
% 
% Anand Sengupta, P. Ajith, 4 February, 2010
% 
% $Id: PsieffDerivative.m 68 2010-01-20 21:58:45Z anand.sengupta $

 setconstants;
 M	= Params.totalMass.*MSOLAR_TIME;
 eta 	= Params.eta;

 switch lambda_i 

   case 'M'
	[dPsiBydM] = dPsiPNbydM(fVec, Params);
	dPsieffBydLamba = dPsiBydM;
	dPsieffBydLamba = M*dPsiBydM;       % returned derivative is w.r.t log(M)

   case 'eta'
	[dPsiBydEta] = dPsiPNbydEta(fVec, Params);
	dPsieffBydLamba = dPsiBydEta;
	dPsieffBydLamba = eta*dPsiBydEta;   % returned derivative is w.r.t log(eta)

   case 'A'
	dPsieffBydLamba = zeros(size(fVec));

   case 'B'
	dPsieffBydLamba = zeros(size(fVec));

   case 'a1'
	dPsieffBydLamba = zeros(size(fVec));

   case 'a2'
	dPsieffBydLamba = zeros(size(fVec));

   case 'a3'
	dPsieffBydLamba = zeros(size(fVec));

   case 'a4'
	dPsieffBydLamba = zeros(size(fVec));

   case 'tc'
	dPsieffBydLamba = 2*pi*fVec;

 end

 return;


%=======================================================================================
function [dPsiM] = dPsiPNbydM(f, Params)
 
 setconstants;
 m	= Params.totalMass.*MSOLAR_TIME;
 eta 	= Params.eta;

 dPsiM = (-37.8569749799331.*f)./power(f.*m,0.6666666666666666) - ...
	55./(384.*f.*power(m,2).*pi) - (27145.*f)./ ...
	(64512.*power(f.*m,1.3333333333333333).*power(pi,0.3333333333333333)) -  ...
	(65.*pi)./(384.*m) + (378515.*f.*power(pi,1.6666666666666667))./ ...
	(96768.*power(f.*m,0.3333333333333333)) -  ...
	(3.183802571673538.*f)./(power(f.*m,0.6666666666666666).*eta) -  ...
	(5.*f)./(128.*power(f.*m,2.6666666666666665).*power(pi,1.6666666666666667).* ...
			eta) - 3715./(32256.*f.*power(m,2).*pi.*eta) -  ...
	(15293365.*f)./(6.5028096e7.*power(f.*m,1.3333333333333333).* ...
			power(pi,0.3333333333333333).*eta) +  ...
	(f.*power(pi,0.3333333333333333))./ ...
	(4.*power(f.*m,1.6666666666666667).*eta) -  ...
	(107.*power(f.*m,0.3333333333333333).*power(pi,0.3333333333333333))./ ...
	(42.*m.*eta) + (38645.*pi)./(32256.*m.*eta) +  ...
	(77096675.*f.*power(pi,1.6666666666666667))./ ...
	(1.6257024e7.*power(f.*m,0.3333333333333333).*eta) -  ...
	(3085.*f.*eta)./ ...
	(9216.*power(f.*m,1.3333333333333333).*power(pi,0.3333333333333333)) +  ...
	(76055.*f.*power(pi,0.3333333333333333).*eta)./ ...
	(221184.*power(f.*m,0.6666666666666666)) -  ...
	(74045.*f.*power(pi,1.6666666666666667).*eta)./ ...
	(48384.*power(f.*m,0.3333333333333333)) -  ...
	(127825.*f.*power(pi,0.3333333333333333).*power(eta,2))./ ...
	(165888.*power(f.*m,0.6666666666666666)) -  ...
	(107.*f.*power(pi,0.3333333333333333).* ...
	 log(power(f.*m,0.3333333333333333).*power(pi,0.3333333333333333)))./ ...
	(42.*power(f.*m,0.6666666666666666).*eta);
 return;


%=======================================================================================
function [dPsiEta] = dPsiPNbydEta(f, Params)

 setconstants;
 m	= Params.totalMass.*MSOLAR_TIME;
 eta 	= Params.eta;
 
 dPsiEta = 3085./(3072.*power(f.*m,0.3333333333333333).*power(pi,0.3333333333333333)) + ...
	(76055.*power(f.*m,0.3333333333333333).*power(pi,0.3333333333333333))./...
	73728. - (74045.*power(f.*m,0.6666666666666666).*...
			power(pi,1.6666666666666667))./32256. + ...
	(9.551407715020614.*power(f.*m,0.3333333333333333))./power(eta,2) - ...
	3./(128.*power(f.*m,1.6666666666666667).*power(pi,1.6666666666666667).*...
			power(eta,2)) - 3715./(32256.*f.*m.*pi.*power(eta,2)) - ...
	15293365./(2.1676032e7.*power(f.*m,0.3333333333333333).*...
			power(pi,0.3333333333333333).*power(eta,2)) + ...
	(3.*power(pi,0.3333333333333333))./...
	(8.*power(f.*m,0.6666666666666666).*power(eta,2)) - ...
	(38645.*pi)./(32256.*power(eta,2)) - ...
	(77096675.*power(f.*m,0.6666666666666666).*power(pi,1.6666666666666667))./...
	(1.0838016e7.*power(eta,2)) - ...
	(127825.*power(f.*m,0.3333333333333333).*power(pi,0.3333333333333333).*...
	 eta)./27648. - (38645.*pi.*log(6))./(21504.*power(eta,2)) + ...
	(107.*power(f.*m,0.3333333333333333).*power(pi,0.3333333333333333).*...
	 log(power(f.*m,0.3333333333333333).*power(pi,0.3333333333333333)))./...
	(14.*power(eta,2)) - (38645.*pi.*...
			log(power(f.*m,0.3333333333333333).*power(pi,0.3333333333333333)))./...
	(10752.*power(eta,2));
 return;


